# Label models (for modelsummary)
modnames = function(l){
  names(l) = paste0("(", 1:length(l), ")")
  return(l)
}

# Function to fix coordinates (Canary Islands)
# Borrowed from https://stackoverflow.com/questions/13757771
fix1 = function(object, params){
  r = params[1]
  scale = params[2]
  shift = params[3:4]
  object = elide(object, rotate = r)
  size = max(apply(bbox(object), 1, diff))/scale
  object = elide(object, scale = size)
  object = elide(object, shift = shift)
  object
}

weekly_pp_int = function(modelnum){

  # Models and checks
  m = lmw_all[[modelnum]]
  ms = lmw_sec[[modelnum]]
  mr = lmw_reg[[modelnum]]
  if(!identical(names(m$coefficients), names(ms$coefficients)) &
    identical(names(m$coefficients), names(mr$coefficients))){
    stop("? -- using lmw_all, lmw_sec, lmw_reg")}

  # Detect org in interaction
  orgname = gsub("week_since_str:", "",
    names(m$coefficients)[grepl(":", names(m$coefficients))])
  print(paste0("Calculating for ", orgname))

  # Main vars
  newdf = expand.grid(week_since_str = seq(1, 15, 1), org = 0:1)
  newdf$patr33_bin = 1
  newdf$sind33_bin = 1
  newdf$CNT_bin = 1
  newdf$UGT_bin = 1
  newdf = newdf[, -which(names(newdf) == orgname)]
  names(newdf)[names(newdf) == "org"] = orgname
  # Rest of variables
  newdf$pop1930_l = mean(data$pop1930_l, na.rm = TRUE)
  newdf$izq1936 = mean(data$izq1936, na.rm = TRUE)
  newdf$analfab30_total = mean(data$analfab30_total, na.rm = TRUE)
  newdf$clerics_all_l = mean(data$clerics_all_l, na.rm = TRUE)
  newdf$prov_name = "toledo"

  ## Predicted probabilities
  # Prepare df
  pp = newdf[, 1:2]
  # All clergy
  pp$y_all = predict(m, newdata = newdf)
  pp$se_all = predict(m, newdata=newdf, se.fit=T)$se.fit
  pp$upr95_all = pp$y_all + qnorm(0.975) * pp$se_all
  pp$upr90_all = pp$y_all + qnorm(0.95) * pp$se_all
  pp$lwr95_all = pp$y_all - qnorm(0.975) * pp$se_all
  pp$lwr90_all = pp$y_all - qnorm(0.95) * pp$se_all
  # Secular clergy
  pp$y_sec = predict(ms, newdata=newdf)
  pp$se_sec = predict(ms, newdata=newdf, se.fit=T)$se.fit
  pp$upr95_sec = pp$y_sec + qnorm(0.975) * pp$se_sec
  pp$upr90_sec = pp$y_sec + qnorm(0.95) * pp$se_sec
  pp$lwr95_sec = pp$y_sec - qnorm(0.975) * pp$se_sec
  pp$lwr90_sec = pp$y_sec - qnorm(0.95) * pp$se_sec
  # Regular clergy
  pp$y_reg = predict(mr, newdata=newdf)
  pp$se_reg = predict(mr, newdata=newdf, se.fit=T)$se.fit
  pp$upr95_reg = pp$y_reg + qnorm(0.975) * pp$se_reg
  pp$upr90_reg = pp$y_reg + qnorm(0.95) * pp$se_reg
  pp$lwr95_reg = pp$y_reg - qnorm(0.975) * pp$se_reg
  pp$lwr90_reg = pp$y_reg - qnorm(0.95) * pp$se_reg

  # Transform
  pp = pp %>%
    pivot_longer(
      cols = c(starts_with("y_"), starts_with("lwr"), starts_with("upr")),
      names_to = c("stat", "clergy"),
      names_pattern = "^(y|lwr90|upr90|lwr95|upr95)_(all|sec|reg)") %>%
    pivot_wider(names_from = "stat", values_from = "value") %>%
    mutate(clergy_lab = recode(clergy, "all" = "All clergy", "sec" = "Secular clergy",
      "reg" = "Regular clergy")) %>%
    mutate(model = orgname)

  # Copy leftist org
  pp$org = pp %>% pull(orgname)

  # Finish
  return(pp)

}

pp_lm_int = function(newdf, model){

  # Get label
  intvar = names(coefficients(model))[grepl("patr33_bin:", names(coefficients(model)))]
  mod_label = gsub("_bin$", "", gsub("patr33_bin:", "", intvar))

  # Prepare df
  pp = newdf[, 1:2]

  # Pred probs
  pp$y = predict(model, newdata=newdf)
  pp$se = predict(model, newdata=newdf, se.fit=T)$se.fit
  pp$upr95 = pp$y + qnorm(0.975) * pp$se
  pp$upr90 = pp$y + qnorm(0.95) * pp$se
  pp$lwr95 = pp$y - qnorm(0.975) * pp$se
  pp$lwr90 = pp$y - qnorm(0.95) * pp$se

  # Mod label
  pp$mod_label = mod_label

  # Copy leftist org
  if("sind33_bin" %in% names(pp)){pp$org = pp$sind33_bin}
  if("CNT_bin" %in% names(pp)){pp$org = pp$CNT_bin}
  if("UGT_bin" %in% names(pp)){pp$org = pp$UGT_bin}

  # Finish
  return(pp)

}

# Add covariates to simulated model
addcovs = function(df, weekly = FALSE){
  df$pop1930_l = mean(data$pop1930_l, na.rm = TRUE)
  df$izq1936 = mean(data$izq1936, na.rm = TRUE)
  if(!weekly){df$dias_control_rep_l = mean(data$dias_control_rep_l, na.rm = TRUE)}
  df$analfab30_total = mean(data$analfab30_total, na.rm = TRUE)
  df$clerics_all_l = mean(data$clerics_all_l, na.rm = TRUE)
  df$prov_name = "toledo"
  if(!"CNT_bin" %in% names(df)){df$CNT_bin = 1}
  if(!"sind33_bin" %in% names(df)){df$sind33_bin = 1}
  if(!"UGT_bin" %in% names(df)){df$UGT_bin = 1}
  return(df)
}

# Function to make First Letter capital
capitalize = function(str){
  c = strsplit(str, " ")
  cu = lapply(c, function(x)
    paste(toupper(substring(x, 1,1)), substring(x, 2), sep="", collapse=" "))
  return(unlist(cu))
}


# custom ggplot theme
my_theme = function() {
  theme_minimal() +#base_family = "Archivo Narrow") +
    theme(panel.grid.minor = element_blank(),
          plot.background = element_rect(fill = "white", color = NA),
          plot.title = element_text(face = "bold"),
          axis.title = element_text(face = "bold"),
          strip.text = element_text(face = "bold", size = rel(0.8), hjust = 0),
          strip.background = element_rect(fill = "grey80", color = NA),
          legend.title = element_text(face = "bold"))
}


# My Stargazer
my_stargazer = function(
  dest_file, model_list, title, label, order,
  omit = "prov",
  covariate.labels,
  notes_table,
  column.sep.width = "-20pt",
  dep.var.labels = NULL,
  dep.var.labels.include = FALSE,
  column.labels = NULL,
  column.separate = NULL,
  se = NULL,
  add.lines=list(c("Province FE", rep("\\multicolumn{1}{c}{Yes}", length(model_list))))){

  filecon = file(dest_file)
  writeLines(
    stargazer(model_list, omit = omit, title = title, label = label,
      order = order, covariate.labels = covariate.labels, notes = notes_table,
      omit.stat = c("f", "ser", "aic", "LL"),
      intercept.bottom = FALSE,
      column.sep.width = column.sep.width,
      multicolumn = FALSE,
      dep.var.caption = "",
      se = se,
      dep.var.labels = dep.var.labels,
      dep.var.labels.include = dep.var.labels.include,
      column.labels = column.labels,
      column.separate = column.separate,
      font.size = "small",
      digits = 3,
      digits.extra = 0,
      star.char = c("+", "*", "**", "***"),
      star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
      notes.align = "c",
      align = TRUE,
      no.space = TRUE,
      add.lines = add.lines,
      notes.label = "",
      notes.append = FALSE),
  filecon)
  close(filecon)

}

# Choose model wisely, instead of [[]]
select_mod = function(modellist, dv, in_f){

  dv_true = unlist(
    lapply(modellist, function(x){
      as.character(formula(x))[2] == dv
    })
  )

  in_f_true = unlist(
    lapply(modellist, function(x){
      grepl(in_f, as.character(formula(x))[3]) 
    })
  )

  select = which(dv_true & in_f_true)
  if(length(select) == 1){
    return(modellist[[select]])
  } else {stop("Problem selecting model, 0 or >1 options")}

}